Source code for hysop.numerics.remesh.remesh

# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.


import sympy as sm

from hysop.constants import __VERBOSE__, __DEBUG__
from hysop.tools.henum import EnumFactory
from hysop.tools.htypes import check_instance
from hysop.numerics.remesh.kernel_generator import Kernel, SymmetricKernelGenerator

Remesh = EnumFactory.create(
    "Remesh",
    [
        "L1_0",
        "L2_1",
        "L2_2",
        "L4_2",
        "L4_4",
        "L6_4",
        "L6_6",
        "L8_4",  # lambda remesh kernels
        "L2_1s",
        "L2_2s",
        "L4_2s",
        "L4_4s",
        "L6_4s",
        "L6_6s",
        "L8_4s",  # splitted lambda remesh kernels
        "Mp4",
        "Mp6",
        "Mp8",  # Mprimes kernels: Mp4 = M'4 = L2_1 and Mp6 = M'6 = L4_2
        "M4",
        "M8",  # M kernels
        "O2",
        "O4",  # Corrected kernels, allow a large CFL number
        "L2",  # Corrected and limited lambda 2
    ],
)


[docs] class RemeshKernelGenerator(SymmetricKernelGenerator):
[docs] def configure(self, n): return super().configure(n=n, H=None)
[docs] class RemeshKernel(Kernel): def __init__( self, moments, regularity, verbose=__DEBUG__, split_polys=False, override_cache=False, ): generator = RemeshKernelGenerator(verbose=verbose) generator.configure(n=moments) kargs = generator.solve( r=regularity, override_cache=override_cache, split_polys=split_polys, no_wrap=True, ) super().__init__(**kargs)
[docs] @staticmethod def from_enum(remesh): check_instance(remesh, Remesh) remesh = str(remesh) assert ( remesh[0] == "L" and (remesh != "L2") or (remesh in ("M4", "M8")) ), "Only lambda remesh kernels are supported." if remesh in ("M4", "M8"): # given M4 or M8 kernels x = sm.Symbol("x") if remesh == "M4": M4 = ( sm.Poly( (1 / sm.Rational(6)) * ((2 - x) ** 3 - 4 * (1 - x) ** 3), x ), sm.Poly((1 / sm.Rational(6)) * ((2 - x) ** 3), x), ) return Kernel( n=2, r=4, deg=3, Ms=2, Mh=None, H=None, remesh=True, P=(M4[1].subs(x, -x), M4[0].subs(x, -x), M4[0], M4[1]), ) else: remesh = remesh[1:] if remesh[-1] == "s": remesh = remesh[:-1] split_polys = True else: split_polys = False remesh = [int(x) for x in remesh.split("_")] assert len(remesh) == 2 assert remesh[0] >= 1 assert remesh[1] >= 0 return RemeshKernel(remesh[0], remesh[1], split_polys=split_polys)
def __str__(self): return f"RemeshKernel(n={self.n}, r={self.r}, split={self.poly_splitted})"
if __name__ == "__main__": import numpy as np from matplotlib import pyplot as plt for i in range(1, 5): p = 2 * i kernels = [] for r in [1, 2, 4, 8]: try: kernel = RemeshKernel(p, r) kernels.append(kernel) except RuntimeError: print(f"Solver failed for p={p} and r={r}.") if len(kernels) == 0: continue k0 = kernels[0] fig = plt.figure() plt.xlabel(r"$x$") plt.ylabel(r"$\Lambda_{" + "{},{}".format(p, "r") + "}$") X = np.linspace(-k0.Ms - 1, +k0.Ms + 1, 1000) s = plt.subplot(1, 1, 1) for i, k in enumerate(kernels): s.plot(X, k(X), label=r"$\Lambda_{" + f"{p},{k.r}" + "}$") s.plot(k0.I, k0.H, "or") axe_scaling = 0.10 ylim = s.get_ylim() Ly = ylim[1] - ylim[0] s.set_ylim(ylim[0] - axe_scaling * Ly, ylim[1] + axe_scaling * Ly) s.legend() plt.show(block=True) # plt.savefig('out/gamma_{}_r'.format(p))